Mon. Not. R. Astron. Soc. 000,[7|fT9l(201 1) Printed 15 November 201 1 (MN KOgX style file v2.2) 



Dusty gas with SPH — I. Algorithm and test suite 

Guillaume Laibe, Daniel J. Price 

Monash Centre for Astrophysics (MoCA) and School of Mathematical Sciences, Monash University, Clayton, Vic 3800, Australia 



15 November 2011 



ABSTRACT 

We present a new algorithm for simulating two-fluid gas and dust mixtures in Smoothed Par- 
ticle Hydrodynamics (SPH), systematically addressing a number of key issues including the 
generalised SPH density estimate in multi-fluid systems, the consistent treatment of variable 
smoothing length terms, finite particle size, time step stability, thermal coupling terms and the 
choice of kernel and smoothing length used in the drag operator. We find that using double- 
hump shaped kernels improves the accuracy of the drag interpolation by a factor of several 
hundred compared to the use of standard SPH bell-shaped kernels, at no additional compu- 
tational expense. In order to benchmark our algorithm, we have developed a comprehensive 
suite of standardised, simple test problems for gas and dust mixtures: dustybox, dustywave, 
dustyshock, dustysedov and dustydisc, the first three of which have known analytic solutions. 

We present the validation of our algorithm against all of these tests. In doing so, we show 
that the spatial resolution criterion A < c s t s is a necessary condition in all gas+dust codes 
that becomes critical at high drag (i.e. small stopping time f s ) in order to correctly predict 
the dynamics. Implicit timestepping and the implementation of realistic astrophysical drag 
regimes are addressed in a companion paper. 

Key words: hydrodynamics — methods: numerical — ISM: dust, extinction — protoplane- 
tary discs — planets and satellites: formation 



1 INTRODUCTION 

Most of our observational information regarding the interstellar 
medium comes to us via dust. Over the last few years, observa- 
tions using the Spitzer and Herschel space telescopes have sub- 
stantially improved our observational sensitivity to and resolution 
of dust emission in a wide range of astrophysical environments. 
Dust grains provide the materials from which the solid cores re- 
quired for the planet formation process are built (see e.g. |Chiang &| 
Youdin 2010a). They also modify the dynamical evolution of the 
surrounding gas by exchanging momentum and energy via micro- 
scopic collisions (Epstein 1924; Baines et al. 1965). Dust grains 
are also the main sources of the opacities in star-forming molecular 
clouds, thus determining their evolution by controlling the thermo- 
dynamics. Accurate determination of both the dynamics of the star 
and planet formation process and its observational signature thus 
require modelling the coupled evolution of gas and dust. 

Given that the N-body evolution of solid particles in a mix- 
ture of gas and solid material would be prohibitive in terms of both 
physical complexity and computational cost, the usual approach 
is to regard the solid phase as a continuum and the mixture as a 
two-fluid system coupled by a drag term. This requires averaging 
physical quantities over a control volume V that is large enough to 
be statistically meaningful but sufficiently small compared to the 
macroscopic scale to allow a continuum description. In diluted as- 
trophysical media, the frequency of collisions between dust parti- 
cles are infrequent enough that the intrinsic pressure of the dust 



phase can be regarded as negligible to a very good level of ap- 
proximation, leaving the dust as a free-streaming collisionless fluid 
whose motion is controlled solely by gravitational forces and the 
drag-term interaction with the gas. 

However, even with a continuous description of the mixture, 
the equations can be solved analytically only for a few simple 
cases (the solutions to two specific problems, dustybox and dusty- 
wave, corresponding to mutually interpenetrating fluids and acous- 
tic waves propagating in a dusty gas, respectively, are derived in 
|Laibe & Price 2 01 la^ As a result, numerical codes have been de- 
veloped in order to model more realistic systems based either on 
A'-body dust particles in Eulerian grid-based hydrodynamics (e.g. 
Fromang & Papaloizou 2006, Paardekooper & Mellema 2006; Jo- 
|hansen et al.|2007||Miniati|2010[|Bai & Stone|2010| > or with a two- 
fluid Smoothed Particle Hydrodynamics (SPH) approach. 

SPH methods for simulating two-fluid mixtures were first de- 
veloped by Monaghan & Kocharyan ( 1995 1, improved (via an im- 
plicit treatment of the drag terms) in Monaghan] {1997} and ap- 
plied in an astrophysical context to the dynamics of dust grains 
in protoplanetary disks (|Maddison et al.||20"03{ |Rice et al.||2004| 
|Barriere-Fouchet et al.|2005 1. The particle-based nature of the SPH 
formalism means that the addition of a dusty fluid is natural. More 
importantly, the drag term that couples the two phases can be im- 
plemented such that the total linear and angular momentum of the 
system are exactly (and simultaneously) conserved, in line with 
the Hamiltonian and exactly conservative nature of the core SPH 
method (e.g. |Price|201 1) . 
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However, the standard methods for treating dusty gas in SPH 
were developed over 15 years ago and our initial attempts to sim- 
ply apply the existing formulations uncovered several issues that 
needed to be addressed. Specifically: 1) the original formulations 
assumed a spatially constant SPH smoothing length; 2) the SPH 
terms for the conservative part of the equations should be derived 
from a Lagrangian; 3) we found that the use of the standard cubic 
spline kernel for drag terms could be significantly inaccurate; 4) we 
encountered several previously unexplored resolution issues in sim- 
ulating two-fluid mixtures; 5) aspects of the implicit timestepping 
scheme suggested by Monaghan ( 1997 1 were found to be problem- 
atic; 6) that treatments of drag have to date generally limited to lin- 
ear drag regimes; and finally 7) that the existing schemes — hav- 
ing been developed with both astrophysical and geophysical dust 
problems in mind — have not been widely benchmarked on prob- 
lems appropriate to astrophysics; Indeed there is a general lack of 
standardised test problems for two-fluid dust/gas codes, a problem 
partially addressed by our first paper (Laibe & Price 201 la). 

In this and a companion paper (Laibe & Price 2011c, here- 
after Paper II), we set out to systematically address issues l)-7) 
in order to develop a robust and accurate code for simulating the 
dynamics of dust in star and planet formation. The importance of 
modelling the dust-gas interaction has been highlighted by recent 
studies showing that instabilities in dust-gas mixtures are good can- 
didates for triggering the concentration of dust during planetesimal 
formation (Goodman & Pindor|2000"l|Youdin & Goodman|2005| >. 

The continuum equations and the relevant parameters describ- 
ing the evolution of dust-gas mixtures are given in Section |2.1| 
Section[2]describes the two-fluid SPH algorithm, addressing issues 
l)-3). The code is benchmarked against a suite of test problems 
that we have specifically designed in order to provide standardised 
benchmarks for other two-fluid gas/dust codes, addressing issues 4) 
and 7) (Sec.|4]l. The implicit timestepping scheme and treatment of 
non-linear drag (issues 5 and 6) are discussed in Paper II. 



2 TWO-FLUID MIXTURES IN SPH 
2.1 Two-fluid gas and dust mixtures 

2.1.1 Densities 

The fact that dust grains of finite size occupy a finite volume is 
accounted for by defining the volume fraction available to the gas 
according to (e.g. Marble 1970; Harlow & Amsde n"["l975) 



1 



Pd 
Pi 



(1) 



This means that the volume densities of gas and dust p e and pd, 
respectively, are distinguished from the intrinsic densities denoted 
p 6 and p d , respectively, according to 



p A = (l-9)pi, 
Pt = (>Pg- 



(2) 
(3) 



The effects associated with finite dust particle size are mostly neg- 
ligible in astrophysical problems since typically the intrinsic dust 
density pd is much higher than the volume density pd and thus 
x 1. We retain these terms, as in earlier SPH formulations (c.f. 
Monaghan & Kocharyan 1995 1 in order to retain a general algo- 
rithm that can be applied both within and outside of astrophysics. 
The conservation of mass in a two-fluid mixture is thus ex- 



pressed by the continuity equations 

djh 
dt 
dpi 
dt 

where v g and \ a are the gas and dust fluid velocities, respectively. 



+ V.(p g Vg) 
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2.7.2 Equations of motion 

The equations of motion, expressing momentum conservation in a 
continuous, inviscid, two-fluid mixture of gas and dust are given by 
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where P 5 and Pd are the intrinsic pressures. Any intrinsic viscosi- 
ties have been neglected. For astrophysical purposes it may be as- 
sumed that the dust is pressureless, i.e. P A = 0. Similarly, the term 
(1 - 6) VP a in the momentum equation for the dust phase — a 
buoyancy term related to the finite size of the dust particles — is in 
general negligibly small. The reader should note that the definitions 
of physical quantities in a two fluid medium require the local fluid 
volume over which the averaging is performed to be defined (see, 
e.g. |Marble|1970)|Fan & Zhu|1998| >. 

The two fluids exchange momentum F dra „, the drag force per 
unit volume, the expression for which is obtained by averaging the 
local drag stress tensor (denoted e^ ag ) over the surface area of the 
dust grains: 



- f e'l dA J . 

yj drag 



(8) 



In the case where the local distribution of dust particles is homo- 
geneous (i.e., dust particles have the same mass, size and intrinsic 
density), Eq.[8]simplifies to 



F drag = *(V g - Yd). 



(9) 



Note that since FY is a force per unit volume, the drag coefficient 
K, has dimensions of mass per unit volume per unit time. This co- 
efficient is related to the drag coefficient on a single grain (denoted 
K s )by 



K 



Pd_ 



a; 



(10) 



where ra d is the mass per grain. The drag force (not per unit volume) 
on a single grain is given by 



"drag 



K s (y g - v d ). 



(11) 



In general K (or equivalently K s ) can itself be a function of the rel- 
ative velocity between the two fluids Av = |v g - Vdl, resulting in 
a non-linear drag regime. In this, Paper I we consider the simplest 
case of linear drag, where K is constant with respect to Av. Exten- 
sion of our scheme to the main non-linear regimes applicable to 
astrophysics are considered in Paper II. 

Finally, it should be noted that in general additional forces 
(e.g. the carried mass, Basset and Saffman forces, Fan & Zhu 1998) 
may be present in two-fluid systems. We have assumed in adopting 
Eqs.|6j|7]that these forces can be neglected for astrophysical appli- 
cations. 
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Figure 1. Computing density in SPH gas (solid points) and dust (hollow circles) mixtures. Standard bell-shaped, Gaussian-like, kernels are adopted (weighting 
indicated by the shading), with a single smoothing length on each particle related to the local number density of particles of the same type. This provides good 
density estimates in both extremes — where dust is concentrated below the gas scale (left panel) and where gas is concentrated below the dust scale (right 
panel). The density of another fluid at the position of a reference fluid (e.g. dust density at the location of a gas particle) is computed using the same smoothing 
length but only neighbours of the desired type. This density is thus allowed to be identically zero, as would be the case for the density of gas-at-dust in the left 
panel (top), or dust-at-gas in the right panel. 



2.1.3 Energy equation 

The evolution equation for the specific internal energy of the gas, 
Kg, is given by 



Pi 



dw g 



-P g [<9V • V g + (1 - 0)V • V d ] + Admg + A.i 



(12) 



where the first term corresponds to the usual compressive (PdV) 
term with the volume reduced by the dust filling factor 6. The sec- 
ond term is the work done by the gas in triggering buoyancy effects. 
The third term is the frictional heating due to the drag force, given 
by 



Adrag = Pg#(Vg - V d ) 



(13) 



The fourth, thermal coupling, term arises when the internal tem- 
perature of the grains differs from the gas temperature (c.f. |Marble| 
1 1970|[Harlow & A msden 1975 ), and in general consists of terms re- 
lated to heat transfer due to conduction (A c011 d) and radiation (A rad ), 
given by 



= A cond + A rad = Q(T g - r d ) + R(aTt 



(14) 



where T g and T d are the temperatures of the gas and dust, respec- 
tively, a is the radiation constant and Q and R are coefficients, de- 
pendent on gas and dust properties, that characterise the heat trans- 
fer. The thermal energy of the dust evolves according to 



Pd ■ 



d«d 
dT 



= -A„ 



(15) 



2.2 Densities for two-fluid mixtures in SPH 

2.2.1 Computing densities in two-fluid SPH 

For two-fluid mixtures, we require a density estimate for each 
phase, corresponding to the exact solution of Eqs. [4] and [5] in SPH. 
The main complication arises from the fact that the local particle 
spacing can be different for each fluid, implying that the two fluids 



should have different resolution lengths calculated based on the lo- 
cal particle number density of their own type. Figure [T] illustrates 
the two limiting cases, i.e. a high concentration of dust in a di- 
luted gas (left panel) and conversely a high concentration of gas in 
a low density fluid of dust (right panel). In each case the smoothing 
length for each type is determined by the local number density of 
particles of the same type. That is, the SPH translation of Eqs. [4] 
and[5]correspond to 



b 

Pi = ^ nijWjjQii); 



h a = n 



hi = n 



l/v 



1/v 



(16) 



(17) 



where v is the number of spatial dimensions and n is a constant 
determining the resolution length as a function of the local particle 
spacing (typically r\ = 1.2 is a good choice for the standard cubic 
spline kernel, see |Price|[20TT) . We adopt the convention that the 
indices a, b, c refer to quantities computed on gas particles while 
(', j, k refer to quantities computed on dust particles. Note that the 
densities and smoothing lengths are independently computed for 
each fluid and are thus — so far — only defined on particles of 
the same type. The numerical solution of Eqs . [T6| and [T7| involves 
determining both p and h for each type simultaneously, since they 
are mutually dependent, thus requiring an iterative procedure. The 
procedure is identical to that adopted in standard variable smooth- 
ing length SPH formulations (see e.g. Price & Monaghan 2007 for 
details). 

An additional complication arises from the need to compute 
the volume filling fraction 6 (Eq. [TJ, defined on a gas particle, a, 
according to 



Pi,a 
Pd 



(18) 



which depends on the density of dust at the gas particle loca- 
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tion. Initially we considered computing this density using a second 
smoothing length for each particle based on neighbours of the other- 
type (this would in turn lead to multiple smoothing lengths on each 
particle if more types were present). However, the key point, illus- 
trated by the right hand panel of Fig.[T] is that it should be possible 
for the local density of dust at the gas location to be identically zero 
(giving 9 = 1) if no dust particles are found within the kernel ra- 
dius computed with the gas smoothing length. Thus, the density of 
dust-at-gas should be calculated according to 



PcU 



" neigh, dust 

^ mjW a j(h a ), 



(19) 



7=1 



where h a is the smoothing length of the gas particle computed us- 
ing gas neighbours as in Eq. [T6] In the case where no dust neigh- 
bours fall within the kernel radius, p ia = 0. This is a very simple 
and efficient method that can easily be generalised to multiple flu- 
ids, requires only one smoothing length per particle and does not 
require any significant additional computational expense. 

The discussion above resolves the first issue highlighted in 
Sec. [T] namely how to deal with variable resolution in multi-fluid 
SPH, generalising the earlier fixed-smoofhing-length formulation 
of Monagh an & Kocharyan| ( |1995| >. A similar discussion to the 
above applies to gravitational force softening on multiple fluids 
in A'-body/SPH codes where the softening formulation is derived 
from a kernel density estimate (Price & Monaghan 20071, in par- 
ticular for the case of a mixture of dark matter and baryonic gas 
(e.g. |Merlin et al.|2010[[Tannuzzi & Dolag|2011J . 



2.2.2 Kernel function 

The kernel function itself can be written as a function of the 
smoothing length h and the dimensionless variable q = |r - t'\/h in 
the form 

o~ 



W(r,h) : 



If 



/(«). 



(20) 



where a is a normalisation constant. The standard Gaussian kernel 
is given by 



/(<?) 



(21) 



where a = tC vI1 . The Gaussian is infinitely smooth (differentiable) 
but has the practical disadvantage of infinite range. A standard al- 
ternative (providing a Gaussian-like kernel but truncated at 2h) is 
the M4 cubic spline kernel (Monaghan 1992 ) 



L 4 (2-qY 

0, 



< q < 1; 

1 < q < 2; 
Q>2, 



(22) 



where a = [2/3, 10/ (In) , in [1,2,3] dimensions. An error 
analysis of the SPH density estimate (e.g. |Price[2 011 ) shows that 
in general the measure of a good density kernel is that the normali- 
sation condition 
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mb w 

W ab 

Pb 



1. 



(23) 



is well satisfied for typical SPH particle distributions, correspond- 
ing to j WdV = 1 in the continuum limit. In general most bell- 
shaped (Gaussian-like) kernels, such as the cubic spline, fulfil this 
criterion ( |Fulk&Quinn|1996) . 

More accurate density estimates can be obtained — at the 
price of additional computational expense — by using kernels with 



extended range that form a better approximation to the Gaussian 
(see Price 2011 1. In particular the M 6 quintic kernel, truncated at 
3h, gives results that are in practice largely indistinguishable from 
the Gaussian, with the functional form 



/(<?) : 



(3 - g) 5 - 6(2 - g) 5 + 15(1 - g) 5 , 0<q<\ 

(3-q) 5 -6(2- qf, \<q<2 

(3 - q)\ 2<q<3 

0, q > 3. 



(24) 



where a = [1/24, 96/(1 199tt), 1 /(20n)]. Use of the quintic is a fac- 
tor of (3/2) 3 as 3.4 times more expensive than the cubic spline (or 
other 2/j-truncated kernels) in three dimensions. 

The functional form of the M4 cubic and M b quintic spline ker- 
nels are shown in the top row of Fig. [2] showing the kernel function 
f(q) (solid/black lines) and its first (dashed/red lines) and second 
(short dashed/green line) derivatives. 



2.3 Equations of motion 

As discussed by Price) |20TT| l, specifying the manner in which the 
density is calculated in SPH can be used to self-consistently deter- 
mine the equations of motion and energy from a variational princi- 
ple, using only the additional constraint of the first law of thermo- 
dynamics. For a two-fluid system, only the dissipationless part of 
the algorithm can be derived in this manner — that is, not including 
the drag terms. 

2.3.1 Lagrangian 

For a system consisting of gas and dust, the Lagrangian is given by 



/> 



1 



-\% - u b (p b , s b ) 



k 



(25) 



(26) 



where u b is the thermal energy per unit mass of the gas (in general 
a function of the entropy i' and intrinsic density p). The equations 
of motion can be derived from the Euler-Lagrange equations, 

d (dL\ _ BL 
dt\dvj ~ ~dr' 

2.3.2 Equations of motion for the gas 

We first consider the evolution of the gas particles. The partial 
derivative of the Lagrangian with respect to the velocity v„ of a 
given gas particle a provides: 



d IdL 
dt\dy a 



' dt 



(27) 



The partial derivative of the Lagrangian with respect to the position 
r a of the gas particle a is given by 



8L 



-z 



in h 



dut 
dpb 



dpb 
dr a ' 



(28) 



where the entropy is constant for a non-dissipative system. Eq.|28| 
differs from the usual expression for single fluids because of the 
distinction between the intrinsic and volume density of the gas 
caused by the finite volume occupied by dust. That is, the ther- 
mal energy depends on the intrinsic density rather than the volume 
density, giving 



du b 
dpb 



Pb 

pi 



Pi ' 



(29) 
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The derivative of the intrinsic density with respect to the particle 
coordinates is given by 



dpb _ dPb 
dr a dp b 

where, fromfT8l we have 



dpb fy>b 
Bb dr a 39 b 



Mb 



(30) 



dpb 
dpt 



Pb 



j_ dpb 
% ®b 98b f 

The spatial derivative of the density sum for the gas (Eq. [T6| is 
given by 

= 7T~ X m f ^ba ~ S ca ) VaWbcQlb), 

dr a a b ^ 
where CI is the usual variable smoothing length term 



(31) 



(32) 



„ _ , dh b v 



dW hc (h h ) 



(33) 



The spatial derivative of the volume filling fraction 9 is given, from 
Eq.^by 

d0 b _ 1 dAu> 
dr a 



Pi dr„ 



(34) 



where 



dpd,b 
dr„ 



■■ ^ mj (d ba - Sja) ^ a Wbj{hb) 



^ m c (S ba - 6 ca ) V„ WbcQib), (35) 



where Q. A is CI computed only using dust particle neighbours, i.e. 



, dhb v-i dW bj (hb) 



dpd,b 



dh h 



(36) 



Collecting Eqs. |28[[36] noting that 6^=0 (since a gas and 
dust index can never refer to the same particle) and using the fact 
&atV e W 6c = -V a W C i, gives 
dL 



dr a 



• - m a ^ m b 

b 

- m a mj 



SlaPl 

P a (i -e a ) 

PaPd.a 



&bPi 



V n W flft Qi b ) 



(37) 



where we have defined 9 to include the correction terms for a vari- 
able smoothing length, i.e. 



9 = 8 + —(1 - 9)(1 - Cl A ). 
Pd 



(38) 



Although this correction is necessary for strict energy conservation, 
it is expected to be negligibly small in practice, since (1 - 9) is neg- 
ligible for small grains and (1 - Cl A ) is 0(h 2 ). Finally, the equations 
of motion for a gas particle, from the Euler-Lagrange equations, are 
given by 

dv n 
dr 



-z 



in,, 



r , l^aWab Ql a ) + ~^V a W ab (h b ) 
"nPS SlbPl 



PaPd.a 



V a W aJ (h a ). 



(39) 



The reader should note that while the first term is a summation over 
gas particle neighbours, the second is summed over dust particle 
neighbours. Eq. [39]may be straightforwardly shown to be a direct 
translation of Eq.[6]into SPH form. Note that the summation over 
dust particles (the buoyancy term) does not involve CI since the 
smoothing length is independent of the dust particle positions. 



2.3.3 Equations of motion for the dust 

The partial derivative of the Lagrangian with respect to the velocity 
v,- of a given dust particle i gives 



d_ldL 
dt\d\i 



dv, 
At' 



(40) 



A buoyancy term arises in the dust because of the dependence 
of the gas internal energy on 9, which in turn depends on the posi- 
tions of dust particles. That is, 

dL du b 
= — > m L 

<5r, 



b d P" 



dpi, 
dr, ' 



(41) 



where 



and in turn, 



dpb 
dr, 



dp b d9 b 
W b dr~' 



(42) 



d9b 
dr, 



1 dpjj, 
Pd 3r, 

— Y^m^u-S^iWbjQib) 

1 ~ D " * Y, m c ( s » - 6 a) Vi W hc (hb) (43) 



Collecting Eqs. |27||43] and noting that 6 b j = 6 ci = 0, we obtain the 
equations of motion for a dust particle in the form 

dv, v P b (l-9 b ), 
dt 



b 



PbPd,b 



-V,Wb,{h b ). 



(44) 



where we have written the kernel using V, W,j, = -V,Wi, to show 
that the force is equal and opposite to that in the gas (Eq.|39|l. It may 
be straightforwardly verified that Eq.|44]is indeed a direct transla- 
tion of Eq.|7]in SPH form. 

Equations [39] and [44] may be combined to show that the total 
momentum is exactly conserved, i.e. 



d 
d7 



0. 



(45) 



2.3.4 Internal energy equation for the gas 

The SPH form of the non-dissipative terms in the internal energy 
equation for the gas (Eq.| 1 2) can similarly be derived from the SPH 
density estimates. In the absence of dissipation the evolution equa- 
tion for a given gas particle a is given by 



(46) 
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Pa dPa 
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dp a 
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dPa 

+ d9 a 
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~d7 ' 


'~ Pi dt'~ 


' Pi 


. dp a 




Pa 


dt 



Using the expressions \3l \ and simplifying using j 1 8ft . we have 

du c 
dt 



9J^dfa + (1 - 9 a )P a dpd,q 

> ; dt p„p A _„ dt 



(47) 



Pa at PaPd.a 

Taking the time derivative of the density sums \16\ and j!9| l we 



have 




dp a 
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dt 


cY a 






dpd,a 
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dt 


i 



(48) 
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(i - n^) 
ci a 



J] m l>(Va-Vb)-VaWab(ha) (49) 
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giving the SPH internal energy equation in the form 



2.4.3 Errors in the integral drag interpolant 



du a G a P a 



At Q a p 2 



2 m b (v„ - v fc ) • V a W ab (h a ) 



(1 ~ 0q)Pq 
PaPi,a 



N neigh, dust 



£ mj(\ a -Vj)-V a W a j(.h a ), (50) 



which indeed can be shown to be an SPH translation of the first two 
terms in Eq,|12| 



2.4 SPH representation of drag terms 

2.4.1 Drag interpolation 
The remaining aspect is to provide an SPH representation of Eq.|9| 



Monaghan & 



specifying the drag term F drag involved in Eqs 
|Kochar yan ( 1995 1 proposed an SPH interpolation of the drag term 
given by 



<A:Av> = v J £(x,x'){[v g (x)- Vd(x')] ■t]tD(x-x',h)Ax', 



where f is the unit vector defined by: 



|x-x'| 



(51) 



(52) 



and v is the number of spatial dimensions of the system (and not the 
inverse of the number of spatial dimensions as one might intuitively 
guess — see below). Monaghan & Kocharyan ( 1995 1 proposed this 
formulation — with velocity difference projected along the line of 
sight joining the particles — mainly because it gives exact conser- 
vation of both linear and angular momentum in the resulting drag 
terms. 

As the SPH interpolation of the drag term does not come from 
the Euler-Lagrange equations derived for non-dissipative term form 
the SPH Lagrangian, the kernel function used in the drag term is not 
constrained to be the same function W used for the density (as as- 
sumed by Monaghan & Kocharyan 1995). Indeed one of our find- 
ings from this paper (discussed below) is that use of a standard 
(bell-shaped) density kernel for drag computations can be signifi- 
cantly inaccurate. We thus use D to denote the kernel employed for 
the drag interpolation. 



2.4.2 Choice of smoothing length in the drag terms 

A key issue is the choice of smoothing length involved in the in- 
terpolation term \5\\ when the gas and dust have different spatial 
resolutions, as illustrated in Fig.[T] We have found from experiment 
that it is very important to smooth the drag term using the maximum 
smoothing length of the two fluids, rather than using an average 
(c.f. Sec. |4.6| and also Ayliffe et al. 2011). Otherwise, unphysical 
resolution-dependent clumping of one fluid below the scale of the 
other can occur. For gas this presents less of a problem because 
there remain pressure gradients that prevent such clumping. How- 
ever, for dust it is crucial since there are no forces that can other- 
wise counterbalance any artificial over-concentration. Since most 
astrophysical problems involve the concentration of dust in a flow 
of gas, a straightforward approach is to simply use the gas smooth- 
ing length when computing the drag interaction. Unless otherwise 
specified (Sec. [431 this is the approach we adopt in this paper. 



The origin of Eq.|5T]can be understood by considering the projec- 
tion of (Av> onto r aa >, the projection of f onto the coordinate a 
(which equivalently denotes the coordinates x, y or z as the system 
is invariant by rotation) and use a Taylor expansion of K and Vd(x') 
around their values on x: 



(KAv)' 



" =v J dx'F 



a:av(x) + 



d(KA\) 



dx 



(x-x') + 0((x-x') 2 )\ 



rD(x-x',h). 



(53) 



giving 



(KAy)" = vKA\(xfl^ + y d^^'V ^ + O (h 2 ) , (54) 



where 



I ap = J Ax'f fDix-x'^h) 
= J dxr^FD(x-x',l 



(55) 
(56) 



This shows that Eq. l |51} is a second-order approximation to the 
drag term, that is, 



(KAv) a = KAV + 0[h 2 ), 
provided the normalisation conditions 

8* 



0. 



(57) 

(58) 
(59) 



hold. Condition {59J and the zeroing of the off-diagonal terms in 
Eq.[58]may be proved straightforwardly by the fact that the integrals 
in i|55|h|56} are odd. The normalisation condition of the diagonal 
terms in Eq.[58]arises because in 3D we have 



/» + py + f~ = 
/' = F> = I . 



1. 



(60) 
(61) 



giving /° = = I~ = l/v. This explains the factor of v in front 
of the drag summation term. 



2.4.4 Discretisation of drag term 

Discretising Eq.[5T|provides the SPH translation of the acceleration 
due to the drag term for both the gas and the dust. Replacing the 
integral by a summation (over particles of the opposing type) and 
pdV with the particle mass, we have 

f ^7 1 = ^ {KAy) = v E m > ^ K • M *«AA»>- (62) 

\ at /drag Pg j PaPj 

for a gas particle and 
(d\i\ 1 



dt 



K 

(KAy) = -v V m b (v bi ■ hi) hiD ib (h b ), (63) 



PbP, 



for a dust particle, where we have defined \ a j = \f, - \j and r a j = 
(r„ - r,)/|r„ - ry|, Importantly, from Eqs. |62(|63| we have: 

( dv„ \ v-i / dv, ^ 



dt 



dt 



0, 



(64) 
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which ensures that the momentum is exactly exchanged between 
the gas and the dust phase by the SPH formalism. Similarly 



I 



m„r„ x 



dv„ 
dt 



+ 



= 0, 



(65) 



' drag i \ / drag 

showing that the total angular momentum is conserved. 

2.4.5 Errors in the SPH drag interpolation 

A key point to note in the formulation of drag terms is that the 
criterion for an accurate kernel drag estimate is different from that 
required for an accurate density estimate (Eq. [23}. Taking Eq. [62] 
and expanding the velocities and the drag coefficient K around the 
position of the gas particle r„, to lowest order, i.e. 



we find 



^1-^ = ^1-^+0(11), 



Pa 



(66) 



(67) 



implying a discrete normalisation condition on the drag kernel of 
the form 

"Hf^/^***- (68) 

, Pi 



The condition {68} implies that the summation on the diagonal 
terms (xx, yy, zz) are equal to unity, while the summations on off- 
diagonal terms (xy, xz, yz) should be zero. The accuracy with which 
this normalisation condition is satisfied depends on the particle ar- 
rangement. While we find that the diagonal terms are well com- 
puted using standard (bell-shaped) kernels, we find that — apart 
from the special case where the dust particles lie on top of the 
gas particles — the off-diagonal terms can be very poorly nor- 
malised. Fig. [3] shows the xy component of Eq. [68] as a function 
of the smoothing length (in units of the particle spacing, Ax) com- 
puted for a dust particle offset by Ax/4 in the x-direction from a 
cubic lattice of gas particles in 3D. Using the cubic spline kernel 
(top left) results in errors of order 5-10% of the diagonal terms for 
reasonable neighbour numbers (h/Ax as 1.1 - 1.5). Furthermore, 
improving the smoothness of the kernel by using the M 6 quintic 
or even the Gaussian (top right) does not significantly reduce the 
error. In numerical tests (Sec. |4.2[ l this manifests as a large error 
in the drag between the two fluids, implying that a more suitable 
kernel is highly desirable. 

2.4.6 Drag kernel function 

After conducting a search for suitable alternative kernels, we found 
that the so-called "double-hump" shaped kernels JFulk & Quinn| 
|1996| > gave a substantial improvement in accuracy — that is, giving 
errors in the computation of Eq. [68] of similar order to the bell- 
shaped kernels in computing Eq. |23| Defining the kernel function 
as previously 



D(r,h)= -g(q), 



(69) 



we construct double-hump kernels from the M 4 cubic and M 6 quin- 
tic kernels using 



giq) = q f(q), 



(70) 



giving, for example, the "double M 4 cubic" (bottom left panel of 
Fig.[2}, the "double M 6 quintic" (bottom right panel of Fig. [2} and 



M 4 cubic spline - 
, , , / 1| — i — . X' . , , , i ," 


" M 6 quintic - 
' . . ■ < rr-r-h-jjii-* . . 1 .' 
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Figure 2. Functional form of the standard bell-shaped cubic spline (top left) 
and quintic (top right) kernels, compared to the double-hump versions of 
these kernels (bottom row). Kernel functions are shown by the solid/black 
lines, while the long-dashed/red and short-dashed/green lines correspond 
to the first and second derivatives, respectively. We find that double-hump 
kernels are significantly more accurate than bell-shaped kernels when com- 
puting SPH drag terms (see Fig.[5J. 



-0.05 
0.1 



\ 1 1 1 1 1 1 

_ \ M 4 cubic spline_ 




i 1 1 

Gaussian_ 








— 1 1 1 1 1— 

double cubic spline 


1 1 1 1 1— r— j 

double Gaussian_ 


1,1,1 


1,1,1, 



1.4 

h/Ax 



1.4 
h/Ax 



Figure 3. Accuracy with which the normalisation condition for the drag 
force is computed using standard bell-shaped kernels (top row) and double- 
hump kernels (bottom row). The plots show the xy component of Eq. |68| 
computed on a dust particle offset from a regular cubic lattice of gas par- 
ticles, as a function of the smoothing length in units of the particle spac- 
ing (h/Ax). With the bell-shaped kernels (top row) the errors are of order 
5 - 10% (the off-diagonal terms should sum to zero). Changing to double- 
hump shaped kernels (bottom row) gives errors < 0.5%. 
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similarly the "double Gaussian". The normalisation constants are 
found in the usual manner by enforcing j DdV = 1, i.e., 



g(q)dV = 1, 



(71) 



where dV corresponds to dq, Inqdq and 4nq 2 dq in one, two and 
three dimensions, respectively. The normalisation constants for the 
double cubic, double quintic and double Gaussian are given by 



C double M 4 
0~ double M 6 
^"double Gaussian 



70 10 

3br' 9n 



1 



42 



1 



60' 277 In' 168/r 



T -v/2 



(72) 
(73) 
(74) 



in [1,2,3] dimensions. The computation of the off-diagonal term in 
( |(iXfr for the double-cubic and double-Gaussian kernels are shown 
in the bottom row of Fig.[3]and indeed show a substantial improve- 
ment, giving errors < 0.5% compared to the 5-10% errors obtained 
using the standard kernels (top row). This improvement in accuracy 
is also reflected in our numerical tests (c.f. Sec. |4.2[ (. 

It is also possible to physically understand the reason why 
the double hump kernel is suited to deal with drag computation. 
In treating multi-fluid interactions, one requires the information of 
one type of particle at the location of a particle of the opposing type. 
Assuming that the number of dimensions of the space is three, the 
SPH smoothing of the physical quantity A corresponds approxi- 
mately to 

^f M q)D( q )d q - 4^rA (g) <?(g+gM) ; <?(g - gM) d g 

Jo Jo * 

_ A (-q M ) + A (q M ) 
2 

(75) 

showing that — for a given particle — the double hump kernel pro- 
vides an average value of a physical quantity stored in the neigh- 
bours of the other species and located at a distance q = q M of 
the particle. For the same reason, the poor accuracy of the bell- 
shaped kernels can be understood because the maximum weight 
corresponds to q = 0, where in general, no particle of the other 
type is present. 



2.4. 7 Frictional heating terms due to drag 

When the system is made of a single gas fluid, the specific thermal 
energy u is a function of state whose total derivative is expressed 
by: 



du = Tds H — -dp. 



(76) 



To generalise this relation with two-fluids interacting with a drag 
term, we derive an additional term for Eq. [76] arising from ex- 
change of momentum (all the other quantities fixed) considering 
a closed thermally isolated system made of gas and dust SPH par- 
ticles, whose energy exchange arises only because of momentum 
exchange (i.e. drag) between two states (denoted i and /, respec- 
tively). Applying the first law of thermodynamics to an infinitesi- 
mal transformation of the system, we have: 



du + de k = (5wj^f + 



(77) 



where u is the total specific internal energy, e k is the macroscopic 
kinetic energy of the system and w is the total work and q is the 



total heat exchanged during the transformation. Assuming that the 
transformation occurs slowly enough for the gas to remain in ther- 
modynamic equilibrium, Eq.[77]reduces to: 



Consequently, 
du| s , p 



du| sp + de k = Tds + -rdp = 0. (78) 
p- 



a a 

_ y (vk + dVkl^J 2 + y v| 
Zj 2 + 2j 2 ' 



As the conservation of the momentum during the transformation 
ensures that: 



we obtain: 



du| s , p = ^ v ka • dv a | SaiP: 



(80) 



(81) 



Using the expression Eq. [62] which gives the evolution of the ve- 
locity of a gas particle due to drag term (i.e. at constant specific 
entropy and density): 



du| s ,p ridu, v 1 



dt 



m-r^r (Vak ■ r a k)r ak D ak (h a ) 

PaPk 



(82) 

which implies that the evolution of the specific internal energy for 
each gas particle is given by: 



/du a \ A drag ^ K ak 

\ dt /dras Pa V PaPk 



Vak ■ hkf D ak (h a ). (83) 



The positive value of d« f ,/dt (and the fact that the kernel function 
is always positive) ensures a positive definite contribution to the 
specific internal energy of each SPH gas particle from frictional 
drag heating. Eq.[83]provides the SPH translation of the drag heat- 
ing term (Eq.|13}. It is straightforward to show that with the above 
expression the total energy is exactly conserved, i.e. 



Zdu a v-i as a v-i 



dv„ 



dv, 
dt 



0. 



(84) 



2.4.8 Thermal coupling terms 



The thermal coupling terms can be expressed in SPH form using 
( |Monaghan & Kocharyan|1995| > 



Atom,* = V m a ^{T a - Tj)W aj + V m a ^a(T 4 a - T*)W a 
p aP j o ' 



for a gas particle, and 



PaPj 



(85) 



A thcl , v = - V m,^(T b - Tj)W bi - V m b -^a(Tt - T?)W b „ 
V P>P b * P»Pi 

(86) 

for a dust particle. Note that we use the standard SPH kernel for the 
thermal coupling terms. Detailed study of the effect of the thermal 
coupling terms, or specific expressions for Q and R, are beyond the 
scope of this paper, and thus for the tests in Sec. [4] we simply set 

T 6 =T.. 
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2.4.9 Drag coefficient 

In general the drag coefficient A" is a function of the properties of 
both the gas and dust. For example in the linear Epstein regime 
relevant to dilute gases in the limit of low Mach numbers, the coef- 
ficient K ak is given by 



Pk Pa 2 

3 \nym A Q a 



(87) 



where s is the grain radius, m& is the grain mass, y is the adia- 
batic index, c s a is the gas sound speed. Since the basic dust-gas 
algorithm described above is insensitive to the specific form of the 
drag, we consider only constant drag coefficients in this paper in 
order to benchmark the method. The detailed implementation of a 
full range of both linear and non-linear physical drag formulations 
is considered in Paper II. 



3 TIMESTEPPING 

3.1 Empirical timestep criterion 

The drag terms impose an additional constraint on the timestep At, 
such that it has to be smaller than a critical value A/ c for an explicit 
scheme (e.g. Leapfrog) to remain stable. Empirically, Mon aghan &| 
|Kocharyan| ( [T995) use the criterion: 



At < min 



(!)■ 



(88) 



which is essentially the minimum of the drag stopping time taken 
over all of the SPH particles. 



3.2 Von Neumann stability analysis 

A more precise criterion can be derived by considering the stability 
of a simple explicit scheme such as the Forward Euler method. We 
consider the evolution of the drag terms over a single drag timetstep 
At, calculating the velocities at the timestep n+ 1 from the velocities 
at the timestep n. Considering only the time-discretisation of the 
equations, we have 



At 



-T-(v g -v d ), 
Pg v ' 



At 



+ — ( V g 
Pd V 



(89) 



(90) 



We then perform a standard Von Neumann analysis, considering a 
perturbation of the velocity field with respect to equilibrium at the 
timestep m corresponding to a monochromatic plane wave, i.e. 



V'V* 



v™ = V™e" 



(91) 
(92) 



where v™ and vjj' are complex constants and k is the wavenumber. 
Substituting Eqs.|91|and|92|into Eqs.|89|and|90|leads to the linear 
system 



v d 



1 - At f At f 
At " 



(93) 



I -At # 

pj 

The two complex eigenvalues A ± _„; of the matrix M are given by 



1 - 







♦') 

Pd/ 


At 




♦*) 


f 


(- 

\pg 


± 2 


(1 


Pd/ 



(94) 



The condition for the numerical scheme to remain stable (|A_| < 1) 
implies a minimum timestep given by 



At < At c 



PgPd 



K(p„ +Pi) 



(95) 



We note that this expression differs slightly to the one suggested by 
Monaghan & Kocharyan ( 1995 1 (Eq,|88[> as it involves the physical 
drag stopping time t s — i.e. the typical time to damp the differential 
velocity between the gas and the dust fluids — rather than - . The 
timestep of Monaghan & Kocharyan ( 19951 is thus correct in the 
limit where the density of one phase is negligible compared to the 
density of the other phase, but becomes erroneous in the case of 
two fluids having densities of the same order of magnitude. Note 
this would apply to grid codes also. 



3.3 SPH explicit timestep 

The stability criterion for the full SPH system (and also for other 
explicit schemes) is expected to be similar to that derived for the 
continuum case (Eq. \95\. The main difference is that the drag co- 
efficient K is in general only defined on particle pairs rather than 
individual particles. We thus take the minimum of Eq. |95|over a 
particle's neighbours, i.e. 



mm 



PaPk 



Ar c ,i 



PbPi 



Kbiipb +pi) 



(96) 



Kakifa + Pk) _ 

for gas and dust particles, respectively. 
3.4 Implicit timestepping 

For strong drag regimes, the timestep restriction imposed by 
Eq. l |96| > becomes prohibitive, and an implicit timestepping algo- 
rithm is required, as proposed by Monaghan (19971. We use only 
explicit timestepping for the tests shown in this paper, with implicit 
timestepping methods discussed in detail in Paper II. 



4 NUMERICAL TESTS 

Despite a number of codes having already been developed for sim- 
ulating astrophysical gas-dust mixtures, none have been bench- 
marked against a wide range of test problems relevant to astro- 
physics. For example, while Monaghan & Kocharyan ( 1995 1, Mad- 
|dison et aL| ( |2003| > consider drag on a single dust particle in a box 
of gas (similar to our dustybox test below), no waves or shocks are 
considered. Similarly Paardekooper & Mellema ( 2006 ) benchmark 
their algorithm against a single dust-gas shock problem with only a 
qualitative solution. Other authors simply check that the timescale 
for settling in an accretion disc is roughly consistent (Barriere- 
Fouchet et al.||2005) or provide no tests at all (|Rice et al.||2004| 



Fromang & Nelson 2005 ; F romang & Papaloizou|2006| >. In the ab- 
sence of known analytic solutions for simple problems, |Johanseri| 
[eTaTI p007) , |Miniati| pOTO} and |Bai & Stone] pOTOl i use the lin- 
ear growth rates for the streaming instability (Youdin & Goodman| 
|2005| > as a test problem, though this is already a complicated prob- 
lem. 

In this paper we present a comprehensive suite of test prob- 
lems designed to investigate all aspects of our algorithm rele- 
vant to astrophysics. These we refer to as dustybox, dustywave, 
dustyshock, dustysedov and dustydisc. Analytic solutions for the 
dustybox and dustywave problems have been derived in Laibe & 
|Price| (201 la) , while the solution for dustyshock is known in the 
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Figure 4. Dust velocity as a function of time in the dustybox problem, using 
2 X 20 3 particles, a dust-to-gas ratio of unity and a constant drag coefficient 
K = 1. Use of the standard cubic spline kernel for the drag terms (short 
dashed/green line) results in errors of order 10% in the velocities compared 
to the exact solution (long dashed/red line). Using a double-hump kernel 
(solid black) improves the accuracy — for the same computational cost — 
by a factor of several hundred, to <, 0.1%. 

limit of high drag. The solutions for dustysedov and dustydisc are 
more qualitative but are important reference problems for astro- 
physical dust-gas mixtures. We consider simulation of the stream- 
ing instability to be of sufficient importance and complexity to be 
covered in detail in a separate paper. Note that all of the tests con- 
sidered in this paper are performed, for simplicity, using a con- 
stant drag coefficient K. Realistic drag regimes are considered in 
Paper II. 

4.1 Code implementation 

Implementation of the algorithm into a standard SPH code with ex- 
plicit timestepping is relatively straightforward. The main changes 
are i) to store a particle type allowing setup and simulation of mul- 
tiple fluids; ii) to compute the densities and smoothing lengths on 
each fluid as described in Sec. |2.2.1| iii) compute the drag term 
between each fluid according to Eqs. |62f|63"1 and the heating term 
given by Eq.|83|and iv) (optional for astrophysics) to implement the 
modifications to the equations of motion due to the volume-filling 
fraction of the dust. We have implemented the two-fluid algorithm 
into both the Af-dimensional ndspmhd test code (Price 2011 1 and 
into the parallel phantom code for 3D problems ( [Price & Federrath] 
[20T0l|Lodato & Price|2010V 

4.2 dustybox: Two fluid drag in a periodic box 

The dustybox problem presented by Laibe & Price| ( [201 la) in- 
volves two fluids in a periodic box moving with a differential ve- 
locity (Avo = v,/,o - v g> o). It is similar to the test performed by |Mon-| 
aghan & Kocharyan ( 1995 1 showing the drag on a single dust grain 
in a box of gas, except that here we consider the dust as a fluid, 



meaning that the densities and smoothing lengths of both phases 
are computed self-consistently. 

4.2.1 dustybox: Setup 

We setup the particles in a 3D periodic domain x,y, z e [0, 1] such 
that the densities p g and p d and the gas pressure P g are constant, 
and neglect the dust intrinsic volume by fixing the volume fraction 
9=1. The box is filled by 20 3 SPH gas particles set up on a regular 
cubic lattice and 20 3 dust particles set up on a cubic lattice shifted 
by half of the lattice step in each direction. The gas sound speed, 
the gas and the dust densities are set to unity in code units and 
no artificial viscosity terms are applied. We give the fluids initial 
velocities v d = 1 and v 5 = 0. 

During the simulation, we verified that both the total linear and 
angular momentum are exactly conserved as expected (Eqs. |64f - 
[65}. We have also verified that 1) the offset of the dust lattice with 
respect to the gas lattice and 2) the timestepping scheme do not 
affect the results. 



4.2.2 dustybox: Choice of drag kernel 

Fig. [4] shows the dust velocity as a function of time in the dusty- 
box test using p g = p d = 1 (i.e., a dust to gas ratio of unity) 
and K = 1, with the exact solution from Laibe & Price (|2011a| > 
shown by the long-dashed/red line. Using the cubic spline M4 ker- 
nel (short-dashed/green line), the errors are of order 10%. Since 
these errors are due to intrinsic bias in the kernel interpolation of 
the drag terms (Fig. [3}, they are independent of resolution, though 
can be improved - at considerable cost - by increasing the ratio 
of smoothing length to particle spacing (i.e., the neighbour num- 
ber). By comparison, use of the double-hump cubic spline kernel 
gives errors < 0. 1% (solid/black line) with no additional overhead 
in terms of cost. 



4.2.3 dustybox: Effect of drag coefficient and dust-to-gas ratio 

Fig.[5]is identical to Fig.|4]but for a range of drag coefficients K = 
0.01,0.1, 1, 10, 100, compared to the exact solution in each case 
given by a solid/black line. Irrespective of the value of K, both gas 
and dust velocities relax to the barycentric velocity (v g = v d = 0.5) 
in a few stopping times 4 = (p g p d )/[A"(p g + p d )]. Using the double- 
hump cubic, an accuracy between 0.1 and 1% is achieved in all 
cases (long dashed/red lines). 

Fig. [6] is similar, but varying the dust-to-gas ratio using 
Pi/pi = 0.01,0.1, 1, 10, 100 (achieved by varying p d with p g = 1) 
and using K = 1 . This changes both the drag stopping time and the 
barycentric velocity towards which the system relaxes. Here again, 
an accuracy between 0.1 and 1% is achieved in all cases. 

4.3 dustywave: Sound waves in a dust-gas mixture 

The exact solution for linear waves propagating in a dust-gas mix- 
ture (dustywave) has been presented by |Laibe & Price| {201 la[ >. 
We have performed a series of tests involving the propagation of 
a sound wave along the x-axis in both one and three dimensions 
in a periodic box, adopting the setup described in Table 2 of Laibe 
& Price ( 201 la I. The dustywave problem is more complex than the 
dustybox problem as the motion of the mixture is driven by both 
the drag and the gas pressure. 
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Figure 5. As in Fig. [4] but using only the double hump cubic kernel 
with a range of drag coefficients K = 0.01,0.1, 1, 10, 100 (top-to-bottom, 
solid/black lines), compared with the exact solution in each case given by 
the long-dashed/red lines. 
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Figure 6. As in Figs. [4] and [5] but varying the dust-to-gas ratio Pd/Ai = 
0.01, 0.1, 1, 10, 100 (top-to-bottom, solid/black lines) and a fixed drag coef- 
ficient K = 1 using the double hump cubic kernel. Exact solutions for each 
case are given by the long-dashed/red lines. 

Specifically, Laibe & Price |2011a| > derive the dispersion rela- 
tion: 

w 3 + iA— + — ) (J - k 2 c 2 co - iK k — = 0, (97) 
\Pg Pi / Pi 

for solutions in the form e '^ kx -" J '\ At high drag, Eq. [97] can be ex- 



panded in a Taylor series, which to first order gives: 

a> = ±kc s -i p& , *^f— ) (98) 

K{p g +p d ) \ 2 ; 

where the effective sound speed is defined according to 

c s = c s A = cJl + yj " . (99) 

The first term of Eq.|98]gives the propagation of the centre of mass 
of the mixture at the effective sound speed c s . The second term 
corresponds to a corrective dissipative term since A e [0, 1], 

4.3.1 dustywave: Setup 

The equilibrium state is characterised by the two phases at rest 
where the gas sound speed, and both gas and dust densities are 
set to unity in code units. In ID this is achieved by placing equally 
spaced particles in the periodic domain x e [0,1]. For the 3D sim- 
ulations, the tests are run in a periodic box x,y, z e [0, 1] with gas 
particles set up on a regular cubic lattice and dust particles set up on 
a cubic lattice shifted by half of the lattice step in each directions. 
As previously, no artificial viscosity is applied. We set the relative 
amplitude of the perturbation to 10~ 4 in both velocity and density in 
order to remain in the linear acoustic regime for which the solution 
in |Laibe &~P rice (201 la) is derived (we have verified that running 
the same simulations setting the relative amplitudes to 10~ 8 gives 
the same results). The density perturbation is applied to the parti- 
cles as described in Appendix B of |Price &" Monaghan (2004). We 
adopt an isothermal equation of state P = c 2 p with c s = 1. 

4.3.2 dustywave: Effect of the smoothing kernel 

The results of the dustywave test in 3D using 2 x 32 3 particles 
with K = 1 and a dust-to-gas ratio of unity is shown in Fig. [7] us- 
ing the standard bell-shaped cubic (green points, lower amplitude) 
and double-hump cubic kernel (black points, correct amplitudes) 
for the drag (left panel) and the double-hump quintic kernel (c.f. 
Sec. |2.4.6| l (right panel), at t = (top) and after 1 and 2 periods 
(middle and bottom panels). The numerical solutions (green and 
black markers) may be compared to the analytic solutions given 
by the solid/red (gas) and long-dashed/red (dust) lines. The ampli- 
tude and frequency of the solution are only correctly captured us- 
ing double-hump kernels. With the double-hump cubic employed 
(left panel, black points) there remains a slight (few %) phase er- 
ror in the numerical solution caused by the remaining kernel bias, 
independent of resolution. A similar error is found generically in 
multidimensional SPH simulations of linear waves (see e.g. Fig. 6 
of Price & Monaghan 2005 ) and in standard SPH is improved by 
using a smoother kernel such as the quintic spline. Indeed, using 
the double-hump version of the quintic kernel for the drag terms 
and the standard quintic for the SPH terms (right panel) we find the 
phase error is smaller by a factor of ~ 5. 

Longer multidimensional simulations of linear waves are 
more complicated in SPH because placement of the gas particles 
on regular lattices are unstable to low-amplitude transverse modes 
that cause the particles to rearrange towards a "glass-like" con- 
figuration (Morris 1996a b). Furthermore, we find that with large 
drag coefficients we require extremely high resolutions to match 
the analytic solutions (see below), which becomes prohibitive in 
3D. In one dimension however, the numerical stability of a sound 
wave is achieved simply by satisfying the courant condition (i.e. 
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Figure 7. Results of the dustywave test in 3D at t = (top row) and after 
1 and 2 wave periods (middle and bottom rows) using 2 X 32 3 particles, 
K = 1 and a dust-to-gas ratio of unity. The analytic solution is given by the 
solid/red (gas) and long-dashed/red (dust) lines. The standard cubic spline 
kernel (green points, left panel) performs poorly on this test. Using the dou- 
ble hump cubic kernel (black points, left panel) both the amplitude and 
frequency are correct but there remains a small phase error due to the ker- 
nel bias which can be corrected by using the smoother double-hump quintic 
kernel (black points, right panel). 



At < 0.3c s /h) and the timestep constraint from the drag (Eq. \95\. 
We thus turn to ID to investigate the full parameter range of the 
dustywave solution. 



4.3.3 dustywave: Resolution requirements at high drag 

Fig.[8]shows the velocity profiles after 10 periods in the ID dusty- 
wave problem for a large drag coefficient (A" = 100) and a dust- to- 
gas ratio of unity, with numerical resolution as indicated. At low 
resolutions (< 256 particles per wavelength) and high drag, the 
amplitude of the wave in the numerical simulations (black solid 
[gas] and open [dust] circles, on top of each other) is severely 
overdamped compared to the analytic solution (red solid and long- 
dashed lines, also on top of each other). This is further illustrated in 
Fig.|9]which shows the kinetic energy as a function of time for sim- 
ulations at different resolutions, compared to the analytic solution 
given by the solid red line. 

Figs. [8] and [9] illustrate a key difficulty that arises when con- 
sidering high drag coefficients, i.e., where the drag stopping time f s 
defined in Eq.|95]is much smaller than the period T of the wave. 
In this case, the drag term efficiently damps the initial differential 
velocity between the gas and the dust in a few t s . However, as the 
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Figure 8. Resolution study for the dustywave test in ID using a high drag 
coefficient (K = 100) and a dust-to-gas ratio of unity using 32, 64, 128, 256, 
512 and 1024 particles from bottom to top. At large drag high resolution is 
required to resolve the small differential motions between the fluids and 
thus prevent over-damping of the numerical solution, corresponding to the 
criterion h < c s t s , here implying > 240 particles. See also Fig. [9] 



pressure continues to drive the propagation of the wave in the gas, 
a small residual de-phasing of order ~ c s t s occurs, which is simply 
the distance travelled by the gas before it is damped by the dust. 
This de-phasing induces a small differential velocity which in turn 
be damped by the drag. This small differential effect dissipates the 
kinetic energy on a timescale ~ t s . 

The spatial de-phasing between the gas and the dust represents 
the smallest length of the problem that must be resolved numeri- 
cally in order to capture the physics of the process. If the spatial 
de-phasing between the gas and the dust is under-resolved, the dif- 
ferential velocity between the gas and the dust is artificially larger 
than the theoretical one, leading to a non-physical over-dissipation 
of the kinetic energy of the system, as observed in Figs.[8]and|9] 

We thus propose a resolution criterion for resolving the differ- 
ential drag of the form 

A < c s t s , (100) 

where A is the resolution length. For SPH, this becomes 

h<c s t s . (101) 

For K = 100 and c s = p„ = p d = 1 in code units this implies 
h < 0.02, i.e. a minimum of ~ 240 particles (assuming rj = 1.2 in 
Eqs.|16|and|17^, which is consistent with Figs.[8]and[9] 

Simulating dust-gas interactions at high drag therefore re- 
quires a high spatial resolution in order to accurately resolve the 
propagation without over-dissipating the energy of the system. This 
can lead to a prohibitive computational cost, somewhat counterin- 
tuitively since the drag simply tends to make the dust stick to the 
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Figure 10. Parameter study of the ID dustywave problem, varying the drag coefficient from 0.01 to 100 (top to bottom) and using two different dust-to-gas 
ratios (Pd/Pg = 1> left figure and Pd/Pg = 0.01, right figure). Results are shown after 5 wave periods using 128 particles except for the K = 100 case where 512 
particles are used to satisfy the resolution criterion h < c s t s . The double hump cubic kernel is used for the drag terms. The results obtained at these resolutions 
are indistinguishable from the analytic solutions (red solid [gas] and dashed [dust] lines) for both the gas (dots) and the dust (circle) particles. 



gas. Most importantly, this requirement is not unique to SPH and 
is a critical issue for any numerical method. Indeed, |Bai & Stone| 
( |2010} find similarly high resolution requirements at short stopping 
times in their simulations of the streaming instability. 



4.3.4 dustywave: Parameter study 

Fig.[l0]shows the results of ID simulations of the dustywave prob- 
lem for 5 drag coefficients (from K = 1CT 2 to K = 10 2 ) and two 
different dust to gas ratio relevant for astrophysical systems (1 and 
0.01 for left and right figures, respectively), showing the velocities 
after five periods compared with the analytic solutions in each case. 
The simulations employ 128 particles except for the K = 100 case 
where 512 particles have been used in order to satisfy the criterion 
{TUT}. For this set of parameters, our method provides results with 
an excellent accuracy (better than one per cent) on the frequencies, 
the amplitudes and the phases of both the gas and the dust velocities 
(and consequently the energy of the system). 

For equal dust to gas ratios (pd/Ps = 1, left figure), both phases 
are equally affected by the drag. At low drag (K = 0.01, top panel 
of left figure), the damping is not efficient enough for gas or the 
dust to be damped as the stopping time is ~ hundreds of periods. At 
intermediate drag (K = 1, middle panel of left figure), the damping 
is the most efficient for the two phases. At large drag regimes (K = 
100), the damping of the differential velocity occurs quickly, but the 
dust density is large enough to distort the gas propagation: the wave 
is de-phased by a half-period compared to the gas-only solution. 

With more typical astrophysical dust to gas ratios (Pd/Pg = 
0.01, right figure), the gas remains essentially unaffected by the 



dust. It thus propagates almost freely in the box at a velocity close 
to the sound speed. By contrast, the dust phase is strongly affected 
by the drag as shown by the K = 0.01 case (top panel, right fig- 
ure), where the damping time for the dust phase is ~ one period. 
The differential velocity between the two phases becomes more 
and more efficiently damped as the drag coefficient increases (right 
panel, from top to bottom), making the dust phase stick to the gas. 



4.4 dustyshock: shock tube in a dust-gas mixture 

Propagation of a shock in a two-fluid dust and gas mixture (the 
dustyshock problem hereafter) has been studied both analytically 
(see e .g. |Rudinger| 1964| > and numerically (see e.g. |Miura & Glass] 
|1982||Saito et al.|2 003 1, using grid based methods. The dustyshock 
occurs in two stages: a transient stage (for which no analytic solu- 
tion is known and therefore studied numerically) followed by a sta- 
tionary stage which consists of the solution for a pure gas solution 
propagating at a modified y and the modified sound speed (Eq.|99| 
see also M iura & Grass|1982] >. In an astrophysical context, simu- 
lations of a dusty shock were used by Paardekooper & Mellema 
(20061 to test their Godunov-type scheme using a Roe Solver de- 
veloped to simulate astrophysical dust and gas mixtures. 

The hypothesis for the dust phase in these seminal studies are 
essentially the same as the ones used in this paper. However, un- 
necessary additional complications arise from their choice of the 
Stokes drag regime (a function of the local Reynolds number for 
the particles), the addition of a heat transfer term (depending on 
the dust conductivity and the Nusselt number of the system) and 
a temperature-dependent gas viscosity. For the purposes of bench- 



14 Laibe & Price 





= 1.5 






-0.5 



Figure 9. As in Fig.[8]but showing the kinetic energy as a function of time in 
the numerical solution at progressively increasing resolution, compared to 
the analytic solution given by the solid black line. The kinetic energy decay 
converges to the analytic solution at ~ 256 - 512 particles per wavelength, 
implying a demanding resolution criterion (h < c s t s ) for high drag. 



marking of our numerical scheme, we instead simulate a simplified 
problem: using a linear drag regime with constant drag term K, no 
heat transfer between the phases and no viscosity other than the 
standard shock-capturing terms used in SPH. While the evolution 
during the transient stage may be different from those considered 
in previous studies, the solution during the stationary stage remains 
unchanged. 



4.4.1 dustyshock: setup 

We setup the dustyshock problem as a two fluid version of the stan- 
dard Sod (19781 problem. Equal mass particles are placed in the 
ID domain x e [-0.5,0.5], where for x < we use p„ = pd = 1, 
v g = Vd = and P a = 1, while for x > p g = p A = 0.125, 
v g = v d = and P g = 0.1. We use an ideal gas equation of state 
P = (y - l)pu with y = 5/3. The density jump means that for SPH 
the resolution is 8 times higher to the left of the shock than to the 
right. We adopt the same initial resolution in both the gas and the 
dust. This differs slightly from the setup used by M iura & Glass| 
( 1982) and |Saito et al.|(2"003[ > where the dust is only placed in the 
right half of the box. Standard artificial viscosity and conductiv- 
ity terms are employed for shock-capturing in SPH as described 
in |Price | ( [20 1 1 [ > with constant coefficients a S p H = 1,/3sph = 2 and 
a u = 1. 



Figure 11. Results of the dustyshock problem with a moderate drag coef- 
ficient K = 1 and a dust-to-gas ratio of unity. This shows the dustyshock 
solution during the transient stage where the analytic solution is not known 
(the solution for the later stationary stage is shown by the dotted red line 
for comparison). Results are similar to those obtained in previous studies 
(Miura & Glass 1982 Saito et al. 2003). Top panels show velocity and 
density in both gas (solid points) and dust (open circles), while bottom pan- 
els show thermal energy and pressure in the gas. Initial particle spacing to 
the left of the shock in both fluids is Ax = 0.001 while to the right it is 
Ax = 0.008, giving 569 equal mass particles in each phase. 



4.4.2 dustyshock: transient evolution 

Fig. [TT] shows the results of a simulation using 2 x 569 particles 
(i.e. a particle spacing of Ax = 0.001 for x < 0) with a moderate 
drag coefficient (K = 1) and a dust-to-gas ratio of unity, show- 
ing velocity and density for both the gas and dust (top panels) and 
the thermal energy and pressure in the gas (lower panels). With 
this choice of drag coefficient the system remains in the transient 
regime at the time shown (since t < t s ). It should be noted that 
while there is no known analytic solution for this stage of the prob- 
lem, the shock profile we obtain is similar to those found previ- 
ously (see e.g. |Miura & Glassl[l982) |Saito et al.|[2003] >. Initially 
as the shock propagates in the mixture, the dust (initially at rest) 
dissipates the momentum and kinetic energy from the gas, lower- 
ing the propagation velocity compared to the ideal (gas only) case 
(dotted red line). The dust density ramps up roughly linearly behind 
the shock, reaching a density near the contact discontinuity roughly 
twice the u nshocked dust density.|Saito et al.|i 



200l) ; |Miura&Glass| 



{1982} and[Paardekooper & Mellema ( 2006) also found a similar 
behaviour, also with a factor of 2 increase in the dust density be- 
hind the shock. The gas-dust interaction to the left of the contact 
discontinuity and in the rarefaction wave has not been previously 
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Figure 12. Results of the dustyshock problem with a high drag coefficient (K = 1000) and a dust-to-gas ratio of unity, thus being in the stationary phase where 
the analytic solution is known (solid/red lines). At low resolution (left figure, same resolution as Fig. 1 1 1 ^ the results are incorrect due to the failure to satisfy 
the resolution criterion at high drag (Eq. |101) . With this criterion satisfied (right panel, using 2x 1 1255 particles) the numerical solution faithfully reproduces 
the analytic result. Thus, as in the dustywave test, extremely high resolution is required to obtain the correct solution at high drag. 



studied since the above authors place the dust only downwards of 
the shock front. We find that the dust density decreases to near zero 
upstream of the contact discontinuity, increasing sharply at the head 
of the rarefaction wave, transitioning smoothly through the rarefac- 
tion wave to match the undisturbed value. We have checked that 
increasing the resolution further does not change the solution sig- 
nificantly for this choice of drag parameters. 



4.4.3 dustyshock: stationary regime 

Fig.|12|shows the results of simulations with a high drag coefficient 
(K = 1000) and a dust-to-gas ratio of unity. In this case, since 
t > t s , the mixture quickly reaches the stationary regime. We are 
thus able to compare the SPH results to the analytic solution given 
by the solid red line (this corresponds to the standard hydrodynamic 
shock solution with modified sound speeds given by Eq. [99}. The 
left figure shows the results at low resolution (2 x 569 particles, as 
in Fig. [TT] while the right figure shows the results at 20x higher 
resolution (2 x 11255 particles). As in the dustywave test, we find 
that at high drag an extremely high resolution is required to obtain 
the correct solution, consistent with our resolution criterion derived 
above (Eq. | 101[ >. If this criterion is not satisfied the numerical shock 
solution is strongly inaccurate (left panel). 



4.5 dustysedov: Sedov blast waves in a dust-gas mixture 

The dustysedov test concerns the propagation of a Sedov blast 
wave in a dust-gas mixture. Although the self-similar Sedov so- 
lution is known for the propagation of a blast wave in a gas phase, 
the solution for a two fluid dust-gas mixture is unknown (though at 
high drag as previously it may be expected that the solution should 
revert to the gas-only solution using the modified sound speed). We 
do not attempt to simulate this problem at high drag as it would 
involve a prohibitive computational expense, instead adopting an 
"astrophysical" dust-to-gas ratio of 1% and a moderate drag coeffi- 
cient such that the presence of dust represents only a small pertur- 
bation to the gas evolution. Results of purely hydrodynamic SPH 
solutions for this test can be found e.g. in|Rosswog & Price ( 2007 ) 
and Springel & Hernquist (20021. 



4.5.1 dustysedov: Setup 

We setup the dustysedov problem in a 3D periodic box (the bound- 
ary conditions are irrelevant for the times shown) at two different 
resolutions, filling the box x,y, z e [-0.5, 0.5] by 50 3 and 100 3 SPH 
particles for both the gas and the dust. Gas particles are set up on 
a regular cubic lattice, with the dust particles also on a cubic lat- 
tice but shifted by half of the lattice step in each direction. We use 
ctsph = 1 and /Jsph = 2 in the artificial viscosity terms, and a n = 1 
in the artificial conductivity term. An ideal gas equation of state 
P = (y - l)pu is adopted with y = 5/3. 
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Figure 13. Cross-section slice showing density in the midplane in the 3D dustysedov problem, for both the gas (left panel) and the dust (right panel) at / = 0.1. 
A dust-to-gas ratio of 0.01 and a drag coefficient of K = 1 have been used with 100 3 SPH particles in each phase. Note the slight difference in the blast radius 
between the dust and the gas, consistent with the response time (f s ) of the dust to the gas drag. 




Figure 14. Results of the 3D dustysedov test, showing the density in the 
gas (left figure) and dust (right figure) from a Sedov blast wave propagating 
in an astrophysical (1% dust-to-gas ratio) mixture of gas and dust with a 
constant drag coefficient K = 1 . The low dust-to-gas ratio means that the 
gas is only weakly affected by the drag from the dust, and is thus close to the 
self-similar Sedov solution (dotted/red line). The dust density is affected by 
the propagation of the blast, resulting in an overdensity that closely mirrors 
the gas overdensity. Results are shown using 50 3 (top panels) and 100 3 
particles (bottom panels). 



In the self-similar Sedov solution, the thermal energy of the 
gas is initially concentrated at r = 0. In the SPH simulation we 
distribute the internal energy of the gas over the particles located 
inside a radius r < r b where r b is set to 2h (i.e., the radius of the 
smoothing kernel). In code units the total blast energy is E = 1, 
with p g = 1 and pj = 0.01. For r > n,, the gas sound speed is set 
to be 2 x 10~ 5 in code units. The dust-to-gas ratio is set to 0.01 to 
be consistent with the value measured for the interstellar medium. 



The drag coefficient is set to K = 1 . Translated to physical units, 
assuming a box size of 1 pc, an ambient sound speed of 2 x 10 4 
cm/s and a gas density of po = 6 x 10~ 23 g/cm 3 the energy of the 
blast is 2 x 10 51 erg and time is measured in units of 100 years, 
roughly corresponding to a supernova blast wave propagating into 
the interstellar medium. Obviously in a real supernova the temper- 
ature inside the blast would be much higher than the sublimation 
temperature of the dust, meaning that it would be quickly evapo- 
rated, so the dustysedov test is mainly useful as a benchmarking 
problem. 



4.5.2 dustysedov; Results 

Fig.[l4]shows the densities of both the gas (left figure) and the dust 
phase (right figure) at t = 0.1 using 50 3 (top) and 100 3 (bottom) 
particles for both fluids. Fig.|13|shows a cross section of the density 
in the midplane in the high resolution (100 3 ) simulation, showing 
the gas (left) and dust (right). As the gas and dust densities are 1 and 
0.01, respectively and the drag coefficient is K = 1 in code units, 
the stopping time is t s = 0.01, which represents 10% of the time 
required for the blast to fill the box. The response of the dust to the 
forcing by the gas drag is therefore of order 10%. Consequently, 
an overdensity in the dust phase forms due to the passage of the 
overdensity in the gas. The overdensities in the gas and the dust 
phases are dephased slightly (seen by comparing the position of the 
peak densities in the gas and dust in Fig. |14) , consistent with the 
finite time (t s ) required for the dust to respond to the gas forcing. 



4.6 dustydisc: Settling and migration in an accretion disc 

Our final test, dustydisc, concerns the evolution of the dust and gas 
mixture in a protoplanetary disc, where vertical settling and radial 
drift of the dust particles (see e.g. |Chiang & Yo udin 20ltib| are 
known to be crucial processes in the early stages of planet forma- 
tion. SPH is well suited to this problem since i) free boundaries are 
trivial to implement and ii) the exact conservation of angular mo- 
mentum by both the gas and dust parts of the algorithm means that 
the problem can be simulated for many dynamical times. We used 
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Figure 15. Rendering of the gas density of a typical T-Tauri Star protoplanetary disc at two different resolutions: 10 5 (left) and 10 6 (right) gas particles. 
Increasing the resolution smoothens the gas phase. The initial dust to gas ratio is 0.01 so that the dust only slightly affects the gas. 



the standard linear Epstein regime given by Eq. [87] The drag force 
is integrated explicitly. The key features — namely both vertical 
settling and radial migration — are expected to occur. We focus 
here on the vertical settling of the grains since the migration is ex- 
tensively discussed in Ayliffe et al. (201 1 ). For this specific test, the 
'artificial viscosity for a disc' described in |Lodato & Price] j2010) 
and implemented in phantom is used. 

n o - 



4.6.1 dustydisc: Setup 

We setup 10 5 gas particles and 10 5 dust particles in a 0.01M Q gas 
disc (with 0.0001M Q of dust) surrounding a 1M Q star. The disc ex- 
tends from 10 to 400 AU. Both gas and dust particles are placed 
using a Monte-Carlo setup such that the surface density profiles of 
both phases are E(r) oc r" 1 . The radial profile of the gas tempera- 
ture is taken to be T (r) oc r~ - 6 with a flaring H/r = 0.05 at 100 
AU. One code unit of time corresponds to 10 3 yrs. A uniform grain 
size of 1 cm is used. 



4.6.2 dustydisc: Resolution issue 

Fig.[T5]compares the evolution of the gas phase, varying the num- 
ber of gas particles from 10 s (left panel) to 10 6 (right panel). The 
number of dust particles does not affect the density profile of the gas 
given the small initial dust to gas ratio. As expected, a smoother gas 
profile is achieved by using a higher resolution. 

Fig.[l6]compares the evolution of the dust phase varying i) the 
smoothing length used to compute the drag term and ii) the number 
of particles in each phase. When the drag term is computed with 
a mixed smoothing length (h = [/?„ + h A ]/2, top left panel), artifi- 
cial structures develop in the dust phase due to over-concentration 
of dust particles below the resolution of the gas. These numerical 
artefacts are removed using instead the gas smoothing length (top 
right panel). Indeed, the gas smoothing length is larger than the 
dust smoothing length since the dust grains concentrate when they 
reach the disc mid plane (see discussion in Sec. |2.4.2) . Smoother 
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Figure 17. Vertical settling of a centimetre dust grain initially located at 
ro = 100 AU and zo = 2 AU (solid/black). SPH results are compared 
to the estimate given by the damped harmonic oscillator approximation 
(pointed/red). The agreement between the numerical and the analytic so- 
lutions indicates that the vertical settling of the dust grain is correctly re- 
produced by the SPH algorithm. The analytic estimate neglects the radial 
drift of the grain and the vertical stratification of the disc. 



dust density profiles are achieved increasing the gas resolution (10 6 
gas particles keeping 10 5 dust particles, bottom left panel). Increas- 
ing the number of gas particles thus reduces the numerical noise in 
the dust phase. Finally, the smoothest dust density profile is natu- 
rally obtained when the resolution in the two phases is the highest 
(2 x 10 6 particles, bottom right panel). 



4.6.3 dustydisc Vertical settling of the particles 

Fig. [TT] compares the vertical settling of a dust particle obtained 
with the SPH simulation and the analytic estimation given by the 
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Figure 16. Rendering of the dust surface density in a typical T-Tauri Star protoplanetary disc with four different configurations: 2 X 10 5 particles using a mixed 
smoothing length (top left) and 2 X 10 5 particles (top right), 10 5 dust and 10 5 dust particles (bottom left) and 2 X 10 s particles (bottom right) using the gas 
smoothing length. Using the mixed smoothing length results in artificial structures in the dust density. Smoother dust profiles are achieved by i) using the gas 
smoothing length and ii) increasing the gas resolution. More accurate results are obtained by increasing the number of both the gas and the dust particles. 



evolution of a damped harmonic oscillator (see e.g. i Garaud & Lin| 
2004). The particle is initially located at r = 100 AU and z = 2 AU, 
and the evolution is computed for 50 code units. The agreement 
between the numerical and the analytic solutions indicates that the 
vertical settling of the dust grain is accurately reproduced by the 
SPH algorithm. It should be noted that the analytic estimation ne- 
glects the radial drift of the grain and the vertical stratification of 
the disc. More precisely, this model assumes an expansion to order 
zero in (-^ /p s )/(ffM/r ) (meaning that the radial and the verti- 
cal motion of the dust particles are decoupled) and to first order in 
Zo/H (the vertical stratification is neglected). The SPH results are 
thus expected to slightly differ from the analytic approximation in 
both amplitude and phase. 



5 DISCUSSION 

In astrophysics, gas and dust mixtures have been predominantly 
studied with grid-based codes. The gas phase is computed as usual 
whereas the dust is treated by using superparticles (e.g. |Youdin &| 
Johansen 2007, Bai & Stone 20101. Computing the drag is usually 
divided into three steps: i) interpolation of the gas velocities at the 
particle positions, ii) calculation of the drag force on the particles 
and iii) attribution of the back-reaction from the particles onto the 
nearby cells. 

Our algorithm has two key advantages compared to the cur- 
rent grid-based algorithms. First, the procedure used for interpo- 
lating the gas velocities at the dust position conserves angular mo- 
mentum exactly, avoiding artificial local torques. Moreover, the in- 
terpolation in current grid-based codes is performed with a stan- 
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dard bell-shaped kernel, regardless of the nature of the error in the 
drag terms. We expect that these aspects of the drag computation in 
grid-based codes may be improved by generalising the techniques 
involved in our SPH algorithm. Second, the equations of motion of 
the mixture (considering the drag terms only) in SPH are invariant 
when permuting the dust and the gas indices, i.e. g «-> d. This sym- 
metry is broken in grid-based schemes where the two phases are 
treated with two different methods (super-particles superimposed 
to a grid). The SPH algorithm provides rigorously identical results 
when interchanging the dust and the gas properties (this has been 
verified on the dustybox problem which involves only the drag, see 
above). 



6 CONCLUSIONS 

We have developed a new general SPH formalism for two-fluid dust 
and gas mixtures, with the aim of simulating the dynamics of dusty 
gas systems in a range of astrophysical contexts. In doing so we 
have generalised the standard methods developed over 15 years ago 
by Monaghan & Kocharyan ( 1995 1 and Monaghan ( 1997) for treat- 
ing dusty gas in SPH. In Sec.[T] we highlighted seven key issues. 
In this, Paper I, we have addressed five of these issues as follows: 
1) we have introduced a simple way to compute SPH densities on 
two fluids with variable smoothing lengths; 2) the conservative part 
of the SPH equations have been derived from a Lagrangian; 3) we 
have demonstrated how the use of "double-hump" shaped kernels 
significantly improve the accuracy of the SPH interpolation of drag 
terms; 4) we find a necessary criterion h < c s t s in order to correctly 
resolve differential motion between gas and dust that becomes crit- 
ical at high drag; we also find it important to ensure /z gas < h iust 
to avoid artificial over-concentration of dust particles, implying a 
higher resolution should be employed in the gas phase relative to 
the dust. 

Finally, to address issue 7), we have presented a compre- 
hensive suite of simple test problems that can be used to bench- 
mark astrophysical dusty gas codes. These consist of the dustybox, 
dustywave, dustyshock, dustysedov and dustydisc problems. The 
first three of these have known (or partially known in the case of 
dustyshock) analytic solutions and can be easily setup in any code 
with standard boundary conditions. We have used these tests to ex- 
plore the issues raised above and have demonstrated that with the 
appropriate resolution criteria satisfied, our formalism is robust and 
provides accurate results. 

The two remaining issues — namely implicit timestepping 
and treatments of astrophysical drag regimes — are addressed in a 
companion paper (Paper II). While this paper concentrates on two- 
fluid gas and dust mixtures, the algorithm is general and can be 
applied easily to the treatment of other multi-fluid systems in SPH 
(e.g. ambipolar diffusion). 
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